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ABSTRACT 

A general three-field variational principle is obtained for the motion of an acoustic fluid enclosed in 
a rigid or flexible container by the method of canonical decomposition applied to a modified form 
of the wave equation in the displacement potential. The general principle is specialized to a mixed 
two-field principle that contains the fluid displacement potential and pressure as independent fields. 
This principle contains a free parameter a. Semidiscrete finite-element equations of motion based on 
this principle are displayed and applied to the transient response and free-vibrations of the couple 
fluid-structure problem. It is shown that a particular setting of a yields a rich set of formulations 
that can be customized to fit physical and computational requirements. The variational principle is 
then extended to handle slosh motions in a uniform gravity field, and used to derived semidiscretc 
equations of motion that account for such effects. 


1. INTRODUCTION 

An elastic container (the structure ) is totally or partly filled with a compressible liquid or gas (the 
fluid}. The fluid structure system is initially in static equilibrium in a steady body force field such 
as gravity or centrifugal forces. We consider small departures from equilibrium that result in forced 
or free vibratory motions. To analyze these motions the fluid is treated as a linear acoustic fluid, 
i.e., compressible but irrotational and inviscid. The purpose of the present work is 

1. To derive variational equations of motion based on a mixed variational principle for the fluid 
subsystem. 

2. To obtain semidiscrete equations of motion following spatial discretization of the coupled prob- 
lem by the finite element method. 

The derivation of the mixed variational principle for the fluid is based on the method of canoni- 
cal equations advocated by Oden and Reddy [13] for mechanical applications. The most general 
dynamical principle derived in this paper contains three primary variables: pressure-momentum 
vector, dilatation-velocity vector, and displacement potential. 
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The general principle is specialized to a two-field functional of Rcissncr type that has pressure 
and displacement potential as primary variables, as well as a free coefficient a that parametrizes 
the application of the divergence theorem. The coupled variational equations are then discretized 
by the finite element method, and semidiscrete equations for a rigid container established. Linkage 
with the structure is then made to establish coupled semidiscrete equations of motion for a flexible 
container. By appropriate selection of the coefficient a a continuum or finite element formulations 
results. One particular setting yields a rich set of symmetric and unsyn.rnetric formulations for 
the transient and free-vibrations elastoacoustic problems. From this set selections can be mac c to 
satisfy various physical and computational criteria. The implications or these selections as regards 
efficiency and numerical stability are discussed. 

The variational formulation is then extented to cover slosh motions in a uniform gravity field 
It is shown that the surface slosh equations may be incorporated as Galcrkin terms in severs 
forms, and that one of these forms merges naturally with the mixed variational principle to form 
an augmented functional. Semidiscretization of this functional produces finite clement equations of 
motions that may be used for a rigid or flexible container. 


2. GOVERNING EQUATIONS 

The three-dimensional volume domain occupied by the fluid is denoted by V. Ihis volume is 
assumed to be simply connected. The fluid boundary 5 consists generally «r two portions 

5 : 5,i U 5,,. 0) 

5 a is the interface with the container at which the normal displacement d u is prescribed (or found as 
part of the coupled fluid-structure problem) whereas 5,, is the “free surface” at which the pressure 
p is prescribed (or found as part of the “fluid slosh” problem). If the fluid is fully enclosed by the 
container, as is necessarily the case for a gas, then 5,, is missing and 5 = 5,j. The domain is rc errec 
to a Cartesian coordinate system (ix, x 2 , X 3 ) grouped in vector x. 

The fluid is under a body force field b which is assumed to be the gradient of a time independent 
potential /?(x), i.e. b = V/3. All displacements are taken to be infinitesimal and thus the fluid 
density p may be taken as invariant. 

We consider three states or configurations: original, from which displacements, pressures and 
forces are measured, current, where the fluid is in dynamic equilibrium at time t, and reference , which 
is obtained in the static equilibrium limit of slow motions. Transient motions are the difference 
between current and reference states. It should be noted that in many situations the original 
configuration is not physically attainable. Table 1 summarizes the notation used in relation to these 

states. 

2.1 Field Equations 

The governing equations of the acoustic fluid are the momentum, state and continuity equations. 
They are stated below for the current configuration, and specialized to the reference configuration 
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Table 1 Notation for Fluid States 


Quantities 


Domain Original Reference Current Transient 


Displacements 

Velocities 

Boundary displacements* 
Displacement potential 
Pressures (+ if compressive) 
Body forces 
Density 


V 

V 

s 

V 

V 

V 

V 


0 

0 

0 

0 

0 

0 

p 


d° 

d° 

d° n 

p° 

b = V(3 


d 

d 

4 

V 

P 

b = V/J 

P 


d = d* -d° 

- . t • 0 

d = d - d 

d n = d^-d% 

\p z= Ip* — ijj 0 

P=P l ~P° 


* Positive along outward normal 


later. The momentum (balance) equation expresses Newton’s second law for a fluid particle: 

pd‘ = — Vp t + b = -Vp‘ + V/?. (2) 

The continuity equation may be combined with the linearized equation of state to produce the 
constitutive equation that expresses the small compressibility of a liquid: 

p‘ = -KVt\ l = -pc 2 Vd‘, ( 3 ) 

where K is the bulk modulus and c = ^/Kfp Hie fluid sound speed. If the fluid is incompressible, 
K,c oo. This relation is also applicable to nonlinear elastic fluids such as gases undergoing 
small excursions from the reference state, if the constitutive equation is linearized there so that 

K = po[dpldp)o. 

The boundary conditions are 

on S,i, p‘ = p‘ <>n S JM ( 4 ) 

where d * is either prescribed or comes from the solution of an auxiliary problem as in fluid-structure 
interaction, and p may be either prescribed or a function of d n and b, as in the surface-wave (“slosh ) 
problem. 

2.2 Integral Abbreviations 

In the sequel the following abbreviations for the volume and surface integrals are used: 

(/)v =' J fdV, Isis =' J s 9dS. M„ 4 = f^sdS, etc. ( 5 ) 
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That is, domain-subscripted parentheses (square brackets) are used to abbreviate volume (surface) 
integrals. Abbreviations for function inncrproducts are illustrated by 


(/.?) 


def 
V = 



(/,,iv« = £ l ! <’ dvdt ' “ L J s j !,dSdt ’ e j' 

3. THE DISPLACEMENT POTENTIAL 


8.1 The Reference State 

Taking the curl of both sides of (2) yields 


curl d =0. 


( 7 ) 


The general integral of this equation for a simply connected domain is 

d‘ = VV>‘ + a + bf, ( 8 ) 

where V’ 1 = 0*(x,t) is the displacement potential, a = a(x) and b = b(x) are time-independent 

vector functions, and t denotes the time. If accelerationless motions (for example, rig. o y 

motions) are precluded by the boundary conditions, a and b vanish. Replacing d = V0 into the 

momentum equation (2) we get _ 

Vp £ = -pVt/>‘ + V/3, 1°) 

which spatially integrated gives _ 

p* = — pi/> 1 + /? + C(t), t ) 

where the scalar < 7(0 is not spatially dependent. Next Integrate the constitutive equation (3) over 

V and apply the divergence theorem to Vd: 

(p‘) v + (pc 2 Vd*)v = (p‘V + [pc 2 d, £ ]s = ^ ^ 


Inserting p‘ from (10) into the above equation furnishes a condition on C{t), which gives 


C(t) = -^WIs + £(*‘)v - \Wv = - vM' + •’*' ~ 


( 12 ) 


££_ 

v ■ V ' ' v v ' ' v 

where o = (l)„ is the fluid volume and 7 = (/)v/« denotes the volume average of a function / 
defined over V. Substituting C(0 into (10) we get 

p* = -p^ 1 -^ 7 ) + (/?-)5) - (13) 

In the static limit of very slow motions, the inertia terms may be neglected and we recover 
reference solution ^2 


p 0 = (/j-/?)-^-K? I*- 




For an incompressible fluid (d n ] s = 0 but c - oo; thus it would be incorrect to conclude that 
pO — p - p. A counterexample to this effect is provided in [14]. 
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3.2 Transient Motions 

Subtracting the constitutive relations at the current and reference states we get 

p = -pc 2 V 2 V> = pc 7 s t ( l5 ) 

where s = -V 2 ^ is called, following Lamb [10], the condensation. Subtracting (H) from (13) yields 

p = -p(0-J)-^Kls. W 

On equating (15) and (16) we get modified forms of the wave equation that account for mean 
boundary surface motions: 


5 = V 2 i/> = + ~[d n ]s, 

c 2 v 


or 


c 2 (v 2 0 - v 2 v>) =-- V> - V’. 


(17) 


The second form follows from 
an incompressible fluid, c — * 
equation V 2 xp = 0. 


-VS = [d n ]s, 
oo and [d„]s 


which is a consequence of the divergence theorem. For 
= 0, and from the first of (17) we recover the Laplace 


S.S Adjusting the Displacement Potential 

If the transient displacement potential is modified by a function of time 


+ P{t), 


(18) 


where tp is the potential of (8)-(17), we may chose P[t) so that c 2 4> - V 2 tf' - s Tor any t. 1 hen 
we obtain the classical wave equation 


c 2 V 2 \p = or 



rp = 0. 


(19) 


In the sequel it is assumed that this adjustment has been made, 
reduces to 

p = -pxp. 


If so, C(l) vanishes and (10) 

( 20 ) 


4. MIXED VARIATIONAL PRINCIPLES 


4.1 Canonical Decomposition 

In this section we derive multifield variational principles for the fluid domain following the canonical 
decomposition method advocated by Oden and Reddy [13]. This method is applicable to self-adjoint 
boundary value problems (BVP) of the form 


Au = / in D 


( 21 ) 
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where u is the unknown function, / the data, A a symmetric linear operator, and D the domain of 
existence of the solution. For time-dependent problems D is the tensor product of the tune domain 
(typically 0 to t) and the volume V. To apply this method, the operator A is factored as 

Au = W*EWu = /, ( 22 ) 

where T and E are linear operators in V and W* is the adjoint of T. This is called a canonical 
decomposition. This decomposition may be represented as the operator composition sequence 

Wu = e, Ee- a, W*a = /, ( 23 ) 

where t and o denote intermediate field variables in D. The three equations (23) are called the kine- 
matic, constitutive and balance equations, respectively, in mechanical applications. The canonical 
representation of boundary conditions on the surface S = 5„ U S„ is 

= g on S u , Bg.o s = h on S„. (21) 

where B s and B* s are surface operators, g and h denote boundary data, and u s = 7 s« 
a s = T s o are extensions of u and a to the boundary S. The extension operators 7s and 6 S often 

involve normal derivatives. 

4.2 The Wave Equation 

The classical wave equation (19) is not a good basis for the canonical decomposition (22). Its 
principal drawback is that the pressure field does not appear naturally as an intermediate variable 
in (23). A more convenient form for our purposes is obtained by taking the Laplacian or both sides 
of (19), and multiplying through by the density p: 

pV 2 (t£ - c 2 V 2 i>) = 0, whence A = pV 2 ^—5 - c 2 V‘^ , / = ( 2; 0 


A suitable canonical decomposition is A — W EW, where 


W = 



E 




W* 




( 20 ) 


in which i = sTl. Boldface symbols are used for W and E because these arc i x l and 4 x 4 
matrices, respectively. The operator product sequence (23) becomes 


e = Wip = 


■ iVip ' 

— 

tv 

s 1 

0 = Ee = 

t/)V0 

-pc 2 V 2 0 

— 

tin 

. V . 


W’ff = pV 2 ^ - pc 2 VU' 


0 . 


(27) 


The intermediate fields e and a are 4 X 1 column vectors. These vectors are partitioned into their 
temporal and spatial derivative subvectors for convenience in subsequent manipulations. Note that 
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the transient pressure p appears naturally as the spatial component of <7. The temporal components 
of e and a are the complex velocity tv and complex specific momentum nn, respective y. 

■ The boundary' portions and S. of (24) are relabeled 5, and 8 pt respectively, to match the 
notation (l). Boundary and initial conditions may be stated as 


= g(x,t) on 5, 1( B*<t(x,() - M x *0 on S i” 

d(x,! 0 ) = d„(x) or m(x,to) = mo(x), d(x,t,) = d,(x) or m(x,f,) = ...,<*)- 


(28) 


Here B end B' are time-independent 4 X 1 and 1 X 4 vectors, respectively related to 
Bs and operators of (24) by B = Bsls and B‘ = B p's, where -» (a scalar) and Ts *4 

matrix) are boundary extension operators for V and <r, respectively. Comparison with (4) 
use of Green’s function reveals that 


B^ = — B5 = [0 0 0 1 


g r = [0 0 0 «f n ], Is = 4^, Ts-l, h ~ P- ( 29 ) 


4.S Three Field Principle 

The most general variational principle for the canonical decomposition ( 2 G) allows the three fields. 
,/, P an d 47 . to be varied independently. The principle may be stated as 6L{rP,e,c) - 0, where 


e, and < 7 , to be varied independe 
functional L is [13] 


L{ 


u,e,e)=L v + ls = t(Ee,e)vx. + (<’. w tt’- v )vxi-(/.t<’)vxi + (»s.Bi/’ gl.Vixi (h.tMs,*^. 


where Ly and Is collect volume and surface terms, respectively. On inserting (27-29) into (3(1) we 
get 

Z.vH( E e, W* — e)vx.= /" f [|e(-v T v + >*■’) - - v) - p( W 4 .)] dV 4f, 

J to J V 

L s ={(Ts,Bxl>-g)s i xt-{h,il>s)s r xi = f \J P(-j£ -dn)dS-i j »’ dL 

J i 0 L J Sj r (3 I ) 


The term (/, VOvxt vanishes and does not contribute to Ly. 

Two Field Principles 

A two field principle of Reissner type can be derived from the functional L by enforcing the iinerse 
constitutive equations e = E' l <7 a prior,’. The resulting principle, winch allows 0 and <r to be 
varied simultaneously, is 6i2(0,<r) = 0, where 

Jt(v, ») = Ry + Rs = - i (E- 0) V X . + (», W *) v X, - (/. *)v X .+ •(».« . D <l> ~ «) s . « < - (* . ** > *, » <■• 
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where Rs — Ls and 


RvW,c) = -i(E-'o,.)vx. + (». W*)y«, = J‘‘ f (^n. T n, - £5 - m T V* - l>V 2 V) W <"■ 

yt0 (33) 

The specific momentum disappears as an independent field if we enforce m = a priori, 

whereupon the functional R becomes a function of rf, and p only and the volume term contracts to 

JM< M = f‘‘ l ( - - "^) J1 ' d ‘- 

To check R = Rvi^iP) + R s we form its first variation* 


(*0 


SR = - + V ! p ,») v>1 - (^3P + + 

^ 3-/. T / \ 1^1 


(35) 


Setting SR = 0 provides the field equations, boundary and initial conditions. 

4-5 Parametrization 

A one parameter family of variational principles can be obtained by transforming all or part of 
the last term in (34), viz. pV 2 0, by the divergence theorem (Green’s first formula for the Laplace 

operator) 


( P Wiv + / WfVpdV- jp^ds- J ' r^dsi + J P^ ds - 

Jy JY J ^ ° * " 


(36) 


Let 0 < a < 1 be the portion of that term to be transformed. Insert pV 2 0 - apV 2 0+ (1 ft)j>V 2 i/> 
in (35) and apply (36) to apV 2 V> to get 


[/ v (ip(V^ )T V^-i^+ a (V^) T Vp-( l-a)pVV) dV-aj^r^dS-af^p d / n dS 


dt. 

(37) 


* The variation of the kinetic energy integral term may be expressed in two different ways, 

J (,v* T ,v#) y „ = (,vV.«^ V)tt - H ('■ v V-.«vv) J„, 

S(pVi T ,V*) vxt = ['H 1 '*],,, • 

depending on whether intention by parts is performed first in lime or .pace, re.pccti.ely. The lirel 
form, which provides physically significant initial conditions, is used in constructing (35). 
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Finally, replace the Laplacian V V left over in (37) by to arrive at the parametrized two-field 

functional! 


JUlM = RaV + Rs = £ l [/ v (- */»(** ) T W - *£5 + a(V^) T Vp - (1 - a )“jf ) ««' 

+ L p[(i-«)^- |i »] <,s+ / s 

S ‘ ' ( 38 ) 

The highest spatial derivative index for both primary variables 0 and p is 1, except if a - 0, in 
which case it is only 0 for p . The two interesting limit cases are of course a = 0 and a = 1 , Tor 

which 

v _ - _ . /•'■ill 

dt, ( 39 ) 


KoW.p) = £ [ j v ' 

Rt (*, P) = 1“ [ X ( - irWfV* - + C ^) T ' Vp) “O' - £ ” d " JS - 1}” - P) dS 


dt. 

(40) 


5. FINITE ELEMENT DISCRETIZATION 


5.1 Discretization of R, t 

In the following we derive semidiscrete finite-element equations of motion based on the R„ functional 
(38). The volume V is subdivided into fluid finite elements. Over each fluid element the state is 
represented by the primary variables \p and p, which are defined as functions of position in the usual 
shape-function interpolation procedure. The finite element interpolation in V may be expressed as 

V>(x,t) = N«,(x) ¥(t), p(x,t) =N p (x)p(t). ( 41 ) 

where 9 and p are computational column vectors that contain nodal values of 0 and p, respectively, 
and N^ and N r are corresponding row-vector arrays of dimensionless shape functions. The specified 
displacement over 5,j is interpolated by 

J„(x,t) = n T d(x,t) = n T N,i(x)d,= NL(x)<1, («) 

where n is the external-normal unit vector on S, h N (l contains the displacement shape functions of 
the enclosing container, N,j„ are these shape functions projected on the outward normal 11 on S . t , 


t If a jt 1, 6R<, = 0 is a restricted variational principle because the substitution V 2 V - <= *0 ho,lls olli y 
at the exact solution. 
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and d contains nodal displacement values. For now the container displacements will be assumed to 
be prescribed, hence the superposed tilde. 

In the following three Sections (5-8) we shall assume that the prescribed-pressure boundary con- 
ditions are exactly satisfied by the finite element interpolation, i.e. p = p on S,,. If so the ,, 
integral of R a simplifies to 

,w) 

which vanishes for a = 1. Inserting (4l}-(42) into the functional (38) with the simplified S v integral 
(43) yields the semidiscrete quadratic form 

i2a (® > p) = -lp* T H®-^p T Gp + a* T Fp + (l-a)[* T Vp-® T Dp f^ T 4]-p T T T d, (44) 


where 


H = [ VNj 

Jv 


VNv, dV = H t , F 


-L 


VNp VNv, dV, 


D = j «- J NjN p dV, V = ^(V. N«) t N p dS, T T = N 


= [ c- 5 NjN, 
Jv 

-I, 


dV = G r , 


In .,„ds, f, 


p V„N,/, dS. 
(45) 

The integration with respect to time is dropped as it has no effect on the variation process described 
below. 


5.2 Continuity Requirements 

The interelement continuity requirements of the shape functions of ip and p depend on the index of 
the highest spatial derivatives that appears in R a . If a ^ 0, this index is 1 for both ip and p an 
consequently C° continuity is required. It is then natural to take the same shape functions for both 

variables: , ^ 

N* = N p ( ,G ) 


with both vectors ¥ and p of equal dimension and evaluated at the same 
matrices in (45) coalesce as 

H = F, G = D = D t . 


nodes. Then some of the 


(47) 


The case a = 0 is exceptional in that no spatial derivatives of p appear. One can then chose C 
(discontinuous) pressure shape functions; for example, constant over each fluid element. If this is 

done, obviously 

N* # Np ( <8 ) 

because ip must be C° continuous. Furthermore the dimensions of p and $ will not be generally 
the same. 
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5.8 Singularity of H 

For later use, we note that matrix H (as well as F if different from H) before the application of atiy 
essential boundary conditions at fluid nodes, is singular because 

He = 0 I 49 ) 

where c denotes the vector of all ones. This follows from (45) and expresses the fact that a constant 
potential generates no pressures or displacements. 

6. TRANSIENT RESPONSE EQUATIONS 


6.1 The Rigid- Container Equations of Motion 

Since R a contains time derivatives of of order up to 2 in *, the appropriate Euler-Lagrangc varia 
tional equation is 


SR a = 


dR a 

3* 


d dR a d‘ 


dRa 


dt d* dr 3 * 


)** + 


dRa 

dp 


6 p-0, 


( 50 ) 


which applied to (44) yields 

[pH* + aFp - (1 - a)Dp + (1 - a)Vp + (1 - 5* = 0 , 

[ - p^Gp + aF T * - (1 - a)D T * + (1 - a)V r * - T T d] 6p = 0 . 

These equations can be presented in partitioned matrix form as 


-(1 - a)D r 0 J 1 P j [ J ~ p * G J l p i l T ^ J 


( 51 ) 


( 52 ) 


where J = (1 — a)V + aF. 

6.8 The Flexible-Container Equations of Motion 

If the fluid is enclosed in a flexible container, the boundary displacements d arc no longer prescribed 
on S d but must be incorporated in the problem by including them on the left hand side of the 
equations of motion. In the sequel, vector d collects all structural node displacements , of which d is 
a subset on S d . Matrix T, suitably expanded with zeros to make it conform to d, becomes T. We 
shall only consider here the case in which the container is modelled as a linear undamped structure 
for which the standard mass /stiffness semidiscrcte equation of motion is 


Md + Kd = f.i + Tp, 


( 53 ) 


where M is the mass matrix, K the tangent stiffness matrix at the reference state, Tp is the 
pressure force on the structure, and f (J is the externally applied force on the structure. Note that K 
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in general mart acecant lor container prestres, effects through the geometric stiffness. Combining 
(52) and (53) we get the coupled system 


M 

0 

0 


* 

' 

0 

pH 

-(l-cr)D 

< 

' 

* 

> + 

0 

-(1 - o)D t 

0 


l P J 

1 


If a = 0, 


'M 

0 

o' 

m 

■ 

0 

pH 

D 

{ » 

► + 

' 

0 

D r 

0 

l p J 



K 

0 

-T 


K 

0 

-T 1 


0 
0 

•T y T 


0 

0 

J T 


-T 

J 

-o~ l G 



d 1 

f ti 1 

\ 

$ 

* 

H 

I 

e-— * 

i 

0 

r 

■ 1 

l p J 

1 0 ) 

( 


-T 

V 



d i 


fu 1 

< 

Of 

1 * 1 

► = < 

-4 


l p J 

1 1 

0 


(55) 


There is little than can be done beyond tins poms, as — — '' T . 

generally different. Although the pressure may be constant over each element, no con ensa ion o 

p is possible in the dynamic case. 

If a = 1, 


'M 

0 



d] 

' 

0 

pH 

0 


• | 

l + 

0 

0 

oj 

1 

[ P J 

1 


K 

0 

-T 1 


0 ' -T 
0 F 
F t -p~ l G 


I 


* > 



(56) 


Note that all these systems, (54) through (56), arc symmetric. 

6.8 Identical Shape Functions 

Further progress in the case a = 1 can be made if we assume, as discussed in §5.2, that the shape 
functions for p and 0 coincide. Taking then (47) into account, (54) simplifies to 


M 

0 

o' 

f d] 

' 

0 

pH 

0 

* 

► + 

0 

0 

0 

l P J 

1 


K 

0 

-T 5 


0 -T 
0 H 
H —p~ l G 


1 r d ) 

til 

< * \ = < 

0 

.Ip) 

1 0 J 


(57) 


The second matrix equation gives p H* + Hp = 0. Since H is nonnegative deli, rite we must have 

P = - P 9. ( 58 ) 

This is the discrete analog of the continuous relation (20) for the dynamic overpressure. For fulutc 
use let us note that if the container is rigid, (57) reduces to 

_p- l Gp + H* =G® +H¥ = T T d. ( 5 «) 

6.4 Unsymmetric Elimination 

If (58) is used to eliminate the pressure vector from (57) we obtain 

SHiHil- 


M pT 
0 G 


{*) 


K 

-T t 


(60) 
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Conversely, eliminating the displacement potential vector gives 


M 

pT T 




(Gl) 


Unlike previous systems, both (60) and (61) are tmsymmeinc. Thus the straightforward elimma 
of a field variable, be it p or t/>, causes symmetry to be lost. These forms will be called unsymmelnc 
two-field forms, or U2 for short. System (60) reduces to (59) if the container is rigid. 

7. REFORMULATIONS OF THE TRANSIENT RESPONSE EQUATIONS 
7.1 SS Forms 

Starting from (57) and (58) it is possible to derive three more symmetric forms that arc formally 
equivalent. One is obtained by differentiating the last matrix equation twice in time, transforming 
the first equation via (57), and finally including (57) premultiplied by p G as third matrix equation. 


M 
P T r 
0 


pT 0 
-pH G 
G 0 



K 0 
0 0 
0 0 


0 

0 

-‘G 



(62) 


0 0 o' 


' d 1 


' -M -pT K 


0 G 0 

{ 

* 

► + 

— pT r pH 0 

4 

0 0 K 

1 

i d J 

1 

-K 0 0 

- 



Another one is obtained by integrating the first matrix equation of (57) twice in time, using (58) to 
eliminate the pressure, and including Kd - Kd = 0 as trivial equation: 


(<W) 


where superposed stars denote integration with respect to f. Finally, differentiating the first matrix 
equation of (63) twice in time, moving pT T d to the left, and including Md < as riua 

equation, we get 


(61) 


The four symmetric forms, (57), (62), (63) and (64), will be called symmetric three field Tonus, 

or S3 forms for short. It should be noted that there is no symmetric S3 form with a state vector 

* * 

consisting of d, p and d. 

7.2 S2 Forms 

Each of the S3 forms has a statically condensable matrix equation that allows ^one field to be 
eliminated. For example, the last matrix equation of (57) is -T d + H® P P w 
can be solved for the pressure vector p if G is nonsingular. Assuming that all matrix inverses 


o o -M 

f d } 

l 

o 

0 

s 

1 

fil 

I 

° 1 

0 pG — pT T 

* + 

0 pH 0 

* 

' " 


-M -pT -K 

l d j 

0 oo 

i <i ) 

1 1 

l -f.i J 
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indicated below exist (more will be said about this later), the condensation process yields four 
two-field symmetric forms: 

pTG-*Hl fell ff,\ 

,hg- 1 hJ\*} W’ 


0 

pH 


[ 


[M 

1° 

M + pTH -l T T 


]{i)*r 


+ pTG _1 T r 
pHG -l T r 


((15) 


GH" 1 T t 


TH 
r‘GH 


[•? :]{?)*[ 


pH + p 2 T T M 1 
pKM _l T 


[ MK"'M „ , 

[pT T K -l M pG + ^T r K- l T 


pMK 


- li 

iTxr-l» 


II). ft ilUl-lAM*-- 


T pT T M 
KM 


(<*) 


(G7) 


(G8) 


These will be called symmetric two-field forms, or S2 forms for brevity. The condensation process 
reduces the number of degrees of freedom but is detrimental to matrix spars, ty. The property 
may be recovered to some exlent by taking advantage of factored forms of the malr.ccs affcclc, > 

the inverses; for example 


— 2 jjsT 


K + pTG 
pHG _1 T r 


pTG _1 H 

pHG" l H 



I T‘ 

K 

0 If K 0 1 

— 

0 H 

0 

pG _1 J [T t II 


(GO) 


Expressions for the matrices in (66)-(68) are given in [3]. 

7 ,s Advantages and Restrictions 

The eight symmetric forms (S3 and S2), plus the two unsymmetric foriris (U2), represent ton 
formulations of the R ,- based fluid-structure interaction problem Tor the .dcnt.cal-shapc-func 
case Although formally equivalent, they may have different behavior in terms of numerica . J 
and computational efficiency. The following items may afTcct the choice among the vanous lorn.* 

Matrix sparseness retention . Matrices G and M are often diagonal. The S2 forms that in ’ 

and If 1 , whether in direct or factored form, arc (other things being equal) preferable o 

Existence of inverses. If the fluid does not have a Tree surface, II is singular on account of (-10), a««I 
consequently (f> 5 )'£es not exist. If the container has some unsuppressed rigid body modes, is 
singular and consequently (€S) does not exist. 

Applied fore, processing. Foots (63) and (67) requite lhat the applied structural forces, f.i, be inte- 
grated twice in time before being used. Doth S2 forms (67) and (68) require add, tonal matrix-vector 
operations on the force vectors. These disadvantages, however, disappear ,n the free-v,bratio„s 

discussed in §8. 

Explicit versus implicit time integration. If M and G are diagonal, both unsymmetric forms (60) 
and (61) are attractive for explicit time integration because the leftmost coefTiccnt matrices 


u 



Table 2 Limit Conditions 


Limit condition 

Matrix 

expression 

Recommended 

form(s) 

Incompressible fluid (c —* oo) 
Cavitating fluid (c -* 0) 

Stiff container 
Hyperlight container 

G -+ 0 
G — ► oo 
K -> oo 

M -* 0 

(60), (61), (62), (66) 
(57), (65) 

(64), (68) 

(64), (68) 


:r^n":;L7wt 

stability or suitability tor perturbation analysts. Of at 1 cond to "» “ contalncd n.,i<l 

ible Ouid case is of central importance. There must b ' » '™ a ' “ C i ' = (6c) <>l)talll lh , 

would behave as a rigid body. Consequently H is nonsmgular. Setting 1 1 

so-called added mass equations (70) 

M u d 4- Kd = l.i i 


where M u is the added mass of the coupled system: 


M a = M + pTH -I T r . 


(71) 


Preservation o, structural rigid tody motions. This k f 

liquid tanks in orbit). covercd in deta il in 58.6. If the fluid is totally 

e p ;Zd by zz:xti:z JL b l ^ *«. w c > *«" ^ -?• 


8 . 


FREE VIBRATIONS 


To obtain the elastoacoustic free-vibrations problem, we make the standard substitutions 

./"* $ = qe JU,t , P = re Jut , f.i = °, 


d = ue 


(72) 


where j = v^I and w is the circular frequency in^the which Tre displayed 

obtain ten algebraic eigenproblems, eight symme . , . thc Appendix. In the following 

below. General properties of these eigensystems are , s “ m ™ eigenvector relations should 

eigenproblem statements, subscript m is a mode index. The following g 


be noted: 


I'm — 


U„,= Wm 5 "™ ( W »» f 


(73) 
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For the unsymmetric forme given in 58.3 one must distinguish between left end 

Supercript R is applied to right eigenvectors wherever necessary; otherw.se left eigenvectors ar. 

assumed. 


8.1 S9 Forms 



(tl. 


0 0 -M 

0 fi G -pT t 
-M -pT -K 



MOO 
0 pH 0 
0 0 0 



(74) 


( 75 ) 


(70) 


( 77 ) 


8.8 S2 Forms 


w; 


M O 
0 pH 


fc)-l K 


+ pTG _1 T r pTG 
pHG -1 T r pHG 


ijJ 


a [M + pTH- l T r Tir'G 
GH _1 T t p _1 GH -1 G 


s _1 nl ( u„, \ 
s-'nj l<J m j’ 

A-.i {:::}• 


*[? k]{--:}=[ pH ;km-^' t «?]{£}• 

*[ 


MK"‘M pMK 

p t t k-‘m P G + p 2 T t 


l T 1 / u m \ [ M Olfu.,,1 

K -, T j ( q,„ J [° ^ II J\‘l»<j 


( 78 ) 

( 79 ) 

( 80 ) ' 
( 81 ) 


8.8 US Forms 


J* 


w- 


M 

pT' 

f Um 1 = 

K 

0 ' 

J l 

0 

G 

1 «im J 

— T r 

H 

l 1m J 


M 0 
pT T G 




( 82 ) 

( 83 ) 
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8.4 Computational Considerations 

The considerations of §7.3 apply for the most part to these ten eigensyslems. However, matrix 
symmetry is more important in free vibrations than in the transient response problem Th.s is 
because eigensolution extraction methods that take advantage or sparsity are more highly developed 
for the symmetric eigenproblem than for its unsymmetric counterpart. For an up-to-date exposition 

of those methods see Parlett [17]. 

The presence of zero eigenfrequencies (w m = 0 roots) may cause serious numerical dim cullies 
in some eigensystem formulations. Two sources of such roots may be distinguished: rigid body 
structural modes, and the constant-potential mode. 


8.5 Rigid-Body Structural Modes 

If the container is not fully supported, Ku r = 0 for structural rigid body eigenmodes u r . If II 
is nonsingular eigensyslems (74)-(75), their condensed versions (78)-(7<>), as well as the two U2 
eigensystems, preserve such modes. To verify this assertion, substitute 


into the Rayleigh quotients (A.12) or (A.15) of the eigensystems. If II is singular, form (79) which 
contains H"\ does not exist, whereas (74) preserves the modes if there exist q r modes such that 
H q r + Tu r = 0. Eigensystems (76)-(77) and (80) do not generally preserve rigid-body modes, 

whereas (81), which contains K , does not exist. 

8.6 Constant Potential Mode and Spectrum Contamination 

Suppose the container is supported so K is nonsingular but the enclosed fluid has no pressure- 
specified surface S P . If so H is singular because of (49). Both U2 eigensystems then possess an 
w = 0 root which conventionally will be assigned modal index 0. This root is associated with the 

following left/right eigenvectors 

Eigensystem (82): Uo = 0, Qo ” e » u o = ^ ^ Cl 

Eigensystem (83): Uo = K l Te, ro = c, u o ~ 

This statement is readily verified by taking the Rayleigh quotients (A.12). The eigenpairs (85-80) are 
collectively called constant potential mode or CPM. The existence and computational implications 
of this mode have been discussed by Geradin et. al. [7]. The mathematical interpretation of (85) is 
“dual” to that of a structural rigid-body mode. Under a rigid-body motion the displacements are 
nonzero but the strains vanish. Under the CPM the potential is nonzero but all fluid displacements 
and dynamic pressures vanish. But unlike rigid-body modes, the CPM has no physical significance. 

it is spurious. 

According to the eigenfunction theory summarized in the Appendix, all non-CPM modes (u,„, 
q mf r m ) of (82) and (83) for m f 0, u m ? 0 satisfy the bi-orthogonality conditions 


<lo = c > 


= e. 


(85) 

( 86 ) 


(0 



= c T (pT T u„. T Gr,„) — C), 


(87) 
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(e T T T K~ l 


M 

P T T ' 

( u "*\ 

0 

G 

\<L»i 


'(T 


t K l Mu, 


+ pT T K” l Tq 


■fit 


+ Gq m ) = 0. (88) 


As regards the symmetric forms, eigensystems 174) ana are - » „ “ ° ... ' 

of H and should not be used. This is because substituting the CP^lcftcgcnvec ,ot W 
one with r m = 0 for (74), produces a Rayleigh quotient for » of the form 0/0. Tins means t hat 
both coefficient matrices have a common null space (the CPM) and every u is an ■ ^ 

a„ eigenproblem h called defective (see Appendix). If one atlempls to „ u mc r, a !y .o v un trc ^ 

defective eigenproblems. nonsensical results can be expected because the whole spectrum likely 

to be contaminated. 

9. SLOSH MOTIONS IN A GRAVITY FIELD 

A liauid with a free surface in equilibrium in a time-independent acceleration field may exhibit sur- 
face waves, informally called -slosh" motions. From an applications standpomt tl.< I most tmpor an 
acceleration fields are gravity and rotational motion, the latter being of ...teres t n Total, ng ^ 
In this section we shall be content with formulating slosh effects in a uniform g V • * 

general fields, including time-dependent body forces, may be variationally treated by . met bo 
of canonical decomposition of the non-homogeneous wave equation, but that general method 
not be followed here as it is not necessary for the gravity case. 

The fluid volume V is in equilibrium in the reference state discussed in §3 1 under the time- 
invariant bod, force per uni, ./ vo/ume b - V0, where 0 is a potential field. As noted above we 
restrict developments here to a gravity field of strength g uniform in space and time. 1 he bo™-!"* 
S is then the egu.TOrium free surface normal to the gravity field. The axes <r„ «, «) selected 
ZlTg acts aLng the -x 3 S -a «is. Hence 0 = -age + B, where B . - -b.trary constant. If 

we chose B so that 0 vanishes at the free surface a - z 0 , then 

0 = -pg{z-z o ). 

In the so-called hydrostatic approximation for small-amplitude gravity waves (9j sloshing is 
ered equivalent to a free surface pressure 

(90) 


where 


J di> - 
r > = dn =fti ° n 


p = p + P9d Ur — P P9 7 ! j 

Here p as before denotes the prescribed part of the pressure (for example, atmospheric pressure) and 
„ is called the elevation of the liquid with respect to the equilibrium free surface This approximation 
assumes that the displacements are infinitesimal and that the 2 -acccleration o t ic s os i mo ion is 

negligible. 

9.1 Variational Principle 

For the variational derivation of “slosh equations" it is advantageous to chose the elevation »; as an 
independently varied field. This choice simplifies the reduction to surface unknowns as well . 
treatment of more complex interface conditions such as capillary effects. 
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To incorporate slosh effects into the mixed variational principles based 1 on the ® j “ 

in §4, it is convenient to follow a Galerkin technique by adding we.glited forms of (88) to 

variation. The following combinations may be considered: 

*(.-#-«..’•£), *(l *{>- 1 sSir.-'’ 1 '),.' 

± (|4-M|^), f . ± (p-p- w .«p) s ,*(I!-". 5 ’)v < 91 > 

V ' / \ (d4> c dip\ 

± (p-p-p»p.«p) Sj ±(^-p.«p) v 

Of these the Urst expression, with signs - and +, oilers two advantage: (1) it is de [’ V * bU [ r °“ 
a functional and (2) it combines naturally with the S„ integral in the first variation (35) 

W parleOised functional the most computationally advantageous choice is again o - 1. 

The expanded functional (40), denoted as R lt , in the sequel, is 

RiniP^’V) = R ' v ~ J tQ [j St P^ dS + f s } P ~^~ P91l) '^ { b>M* dS \ dL (92) 

where R v v is the volume integral of (40). Note that setting q - 0 restores R t . 

9.2 Finite Element Discretization 

In addition to the assumptions (41), (42) and (46) we interpolate i, as 

rj = N„q on 5,„ 

where column vector , contains „„ field elevations at nodes on S and row vector N 
corresponding elevation shape functions. The semidiscrete quadratic form for (92), again excluding 

the time integral, is 

*.,<*,*•) - - i ;P T <W(H - Q„ +) * -P^ T a + - *■»> - ^ 

where 

q„ + = £ NjVN*dS, Q„ + = ^NjVNodS, S = J $ NjN , dS = S T , f* = /^ VN Jp. 

matrices are partitioned as 

. i r ty t t 1 

, (9G) 


(03) 


*={*;}’ p= {pi}’ q - + - |q " 01 • Q,+ 


Q„ o 
o o 


ii = 


ii„ II.. 

II.. H,m, 
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whore *. contains potentials at »„* nodes of elements connected to S„ and p. contains 
Is. The dimensions of Q„ and Q„ are n„ x In general „„ < »„u , (m « *. *» d cl » 
Also typically n„ « n* = n„ as the latter pertain to a volume mesh. If u » tnte, pointed by 

• same surface functions as p, i.e. N,, = N,, on S,„ then 


Q n = Q„ = Q. Qr,+ = [ Q °1» Q 


_ q 0 
1*+ “ [ 0 0 


( 97 ) 


9.3 The Rigid Container 

The following equation, of motion for the rigid but mobile container are obtained on rendering (91) 
stationary: 

* 

P 

Pff*? 


pH 

0 

o’ 

f • 1 

0 

0 

0 

{ p + 

0 

0 

0 

v p#*i J L 


0 

H - Q,,+ 
Q„ + 


H - Qp+ 
-p-‘G 




0 


0 

-s 



V. 


( 98 ) 


matrix equation _ _ _ ^ (99) 


pH# + (P + R+)® = 4 + P( 11 ~ Qr+) G lT l1. 


where 


iT c — i-t 


R + = Q^ + S 


ir,+ 



[ q t s _1 q ol 


'© 

Pi 


r 

O 

o 


r 

o 

o 


= rJ, p = p(h-qJ + )g- 1 (ii-Q„ + ) = p ■ 

( 100 ) 

, _ . , a i ..r c that is n For most real liquids, acoustic and 

rrtii^tvr 

™"^bo™ ” II the dhpiacemcut-potentia, subspace dor, nod 

by the second matrix equation of (98): 

(H-Q,, + )» = Td. (l01) 

For simplicity let u, assume that the container is not only rigid hut motionless, that is, d = 0. The 
incompressible-fluid equations become 


hi 


H„ 

H„„ 


{*:)*[! :]!::!■{■)■ 


(102) 


. t • a f Ti n 1* - 0 Subvector 9 V may be statically condensed from these 

subject to the constraint (H - Q,,+)* - u - auovecior / 

two relations, which may be combined as the system 

[*• o]{!;} + [h. r q. H, o Qr ] {*;}“{ «}■ . (,03) 
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where A* are Lagrangian multipliers (in fact, the pressures at nodes of *.), »"<> 

H. = H„ Q,= [o]' 

If d # 0 the force term in (103) must be appropriately modified. 

Q.j The Flexible Container 


(104) 
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0 

0 

0- 


f n [ 

0 

pH 

0 

0 

< 

* + 

0 

0 

0 

0 


p 1 

. 0 

0 

0 

0. 


\P99 ) 


K 

0 

— T t 
0 


0 

0 

H - Q„+ 

Q*+ 


-T 

H-C&. 
— P -1 G 
0 


0 

0 

-s 



d 1 


rf.n 


* 

I P 1 

i = 4 

U 1 

Ini 


^ pgn J 

1 

l 0 ) 


(105) 


Eliminating t? and p by static condensation yields 


M 

0 


0 

pH 


-Y 
P + 


yiii-ui- 


where 


(106) 

(107) 


K = K + pTG 1 T t , Y = pTG l (II Qj>+)- 

System (106) is the counterpart of (65). If the fluid is treated as ineompressihlc, a subspaee reduction 
procedure similar to that used in §9.3 can be invoked. 

9.5 Slosh Vibrations , ... 

a*- 

Z 'ML compressible field held in 

eigenproblem associated with (103), suppressing the modal index m for s.mphcty, m y 

(108) 


as 


u 


pH, 

0 




L 1 . V, A „ „r © and A h respectively. The last matrix equation 

ravSK* " - b r iz z amplit " dC! ° r * ’ ‘ 

r, = Using these relations we can transform the e lg enproblem (108) to 

UHcA qT o c ]{^ 


P9 S 
0 


i 


in which C = QH7 l Q T (ll0) 

and r. are Lagrange-multip.ier modal amplitudes 

;Cw™d1o°rt: wE* lin^ase. This Knife element reduction-Uvsnrface technique provides 
an alternative to boundary integral methods (1,8). 
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10. CONCLUDING REMARKS 


Displacement-potential formulation, are of practical interest in lloid-struclure transient-response 
uispiacement, v * effective numerical computations. For some 

and vibration an “ y ^“ * T j 2 14 ] d rerc , ences therein. The preceding treatment unifies a 

ZSZtZZZ emttinuuns-based and aigebraic statements “'"> W 

problem. It may be further extended in the following directions. 

(11 The inhomogeneous wave equation .’V 1 * -* = /./ * l,,e bodj ' r °; ce “ 

’ i, time-dependent and V2b / 0. Additional forcing term, appear in the equations of motion. 

These are of interest for slosh of fluids in rotating containers. 


(2) Retaining the specific momentum m as independent field in functional (33) 

(3) Inclusion of additional physical effects: capillarity, cavitation and viscosity. 


Acknowledgements 

The work of the first author was supported by NASA Lewi. Research Center under Grant NAG MM 
monitored by Dr. C. C. Cham». 


REFERENCES 

1. DeRuntz, J. A. and Geers, T. L., Added mass computation by the boundary integral method, hit. J. 
Hum. Mtih. Engrg., 12, 1978, pp. 531-550 

2 Felippa, C. A. and DeRuntz, J. A., Finite element analysis of shock induced hull cavitation, Ompuirr 
Methods in Applied Mechanics & Engineering, 44, 1984, pp. 297 337 
3. Felippa, C. A., Syuim.tri.ation of the contained compressible laid vibration cigenprobleni, Comm. Appl. 
burner. 1, 1985, pp» 241-247 

5. Felippa, C. A., Symmetrisatlon of coupled sigenproblem. by eigenvector align, citation, Comm. A, fL Nu- 

mer. Meth., 4, 1988, pp. 561-563 

_ T a t? ,vv, r C Finite-element /boundary-element analysis of multiple structures 

5 f;er,:kl,tr.Lr:^iic wat F ,:Nut, i ca| / Mo,hod, for LVansicd and Conpied ed. 

by R W. Lewis e t ol., Pmeridge Press, Swansea, U. K., 1984, pp. 

7 Geradin, M., Roberts, G. and Hack, A., Eigenvalue analysis and transient response of 11., id structure 
interaction problems, Eng. Comp., 1, 1984, pp. 151-160 

8. Khabbaz, G. R., Dynamic behavior of liquid in elastic tanks, AIAA Journal, 0, 1970, pp. 1985-1990 

9. Kinsman, D-, Water Wa^s, Dover, New York, 1965 

10. Lamb, H., Hydrodynamics, Dover, New York, Gth ed. 1945 

r " • 


22 


ASME/AMD Volume 64, 1984 

,3. • Oden , J. T. and Redd,, J. N., V.ri.tional Method, in Tl, eorelie.l Mcel, anion, Snd cd., Springer- Verlag, 
Berlin, 1983 

H - 2S LTaL v r^ cd - by 

R W Lewia, P. Bettess and E. Hinton, Wjley, 1984 

10. Oh., on, R-, .nd Felippa, C. A., The effect ol wall motion on the governing eqnat.on. oi contained 
fluids, submitted to J . Applied Mechanics 

17. Parlett, B., The Symmetric Eigenvalue Problem, Prentice-IIall, Englewood ClilTa, N. J., 1980 


23 


Appendix A: TIIE GENERALIZED ALGEDRA1C EIGEN PROJJLEM 


Some facts about the algebraic eigenproblem are collected here for convenient reference. These facts are 
relevant to the study of the free vibrations of the coupled fluid-structure system. 

A.l The Standard Unsymmetric Eigenproblem 

The standard eigenproblem for a real unsymmetric square matrix A may be stated as 


Ax* — A,x,-, 


(A.l) 


where A, the eigenvalues (which may be complex), and x, the corresponding right eigenvectors normalized 
to unit length. The eigenproblem for the transposed matrix is 

A T y i = A.y,-. ( A - 2 ) 

This problem has the same eigenvalues but in general the eigenvectors y,- will be different. The y, are called 
left eigenvectors of A because they satisfy the problem yfA = A,y t ; this in turn explains the qualifier 
‘right’ applied to x,-. The system of left and right eigenvectors of A satisfies bi-orthogonality relations: 


T 

y i xy 


r o if * ^ j, 
l P. ift=j. 


(A.3) 


This fit is called the condition number of A, with Tespect to the eigenproblem (A.l); it is always less or equal 
than 1 in absolute value, and may be zero in pathological cases. (The closer to 1, the better conditioned A, 
is.) 

Premultiplying (A.l) by y, and assuming that ^ 0 yields 


A *• =yfA Xi/p,- = xj A T yJm, (^-'0 

which is the Rayleigh quotient for unsymmetric matrices. If /*■ = 0 and y] Ax = 0, (A. 5) takes the 
undetermined form 0/0 so every A,- is an eigenvalue. In such a case the eigenproblem (A.l) is said to he 
defective. 

A.£ The Standard Symmetric Problem 

If A is symmetric then x, = y;» m = 1 and (A.3) reduce to the usual orthogonality conditions 


T 

X, X 




if * ? 
if » = j. 


(A.5) 


whereas (A.4) becomes the usual Rayleigh quotient for a unit length vector: 

A, = xf Ax*. 


(A.6) 
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A. 3 The Generalized Unsymmelric Etgenproblem 
The generalized unsymmetric eigenproblem is 

Ax, = A.Bx,, 


(A.7) 


where A end B ere ..^metric reel melricee. Aeseminj the! B“ exiete. thie problem een be redeeed to 
the standard problem 

Cx,- = A.Xi, 


(A. 8) 
(A. 9) 


in which C = B" l A. The transposed problem is 

C T », = A T B -T i; = A i.i. 

Defining B T y, = *«> ( A - 9 ) can be tran9fomed to 

A T y, = A.B T y t . 

The bi-orthogonality conditions (A.3) become 

_ _ - T ( 0 iti yt j, 

x j = y, — x ' ® — \ m if t = j. 

The Rayleigh quotient (A.4) generalizes to 

_ yf Ax - _ y^ Ax » 

1 " yTBx, m 

A in KA 1 if (A 121 takes on the form 0/0 for some t, every A,- is an eigenvalue and the eigenproblem ( A - 7 ) 

a»=2= 

the entire spectrum. 

A-4 The Generalized Symmetric Eigenproblem 

If both A and B are symmetric, (A 131 

x,=y,, *• = B y,- 

and we recover the usual orthonormality conditions 


( A..10) 


(AH) 


(A. 12) 


T-rj f 9 it i j, 
C ’ Bx ' = U if » = 3- 


(A.I4) 


In mechanical vibration problems for which B is the mass matrix, /i; is called the generalized mass. Finally, 
(A. 12) reduces to the usual Rayleigh quotient 


xf Ax, 

’ " xTBx. 


(A. 15) 
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